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Abstract 

In this work a novel method for the analysis with trimmed CAD surfaces is presented. 
The method involves an additional mapping step and the attraction stems from its sim¬ 
plicity and ease of implementation into existing Finite Element (FEM) or Boundary 
Element (BEM) software. The method is first verified with classical test examples 
in structural mechanics. Then two practical applications are presented one using the 
EEM, the other the BEM, that show the applicability of the method. 
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1. Introduction 

The original idea of isogeometric analysis was to use data from CAD programs 
directly for describing the geometry of the problem without the need to generate a 
mesh. Geometry data are provided by CAD programs in a standard ASCII format 
(IGES standard^) and can be read by the analysis program. CAD programs describe 
surfaces using NURBS functions. Eor complex surfaces, especially involving a union 
of surfaces, trimming is used. The data provided in the IGES file contain NURBS 
information of surfaces (control points, knot vectors and weights), as well as trimming 
information, if the surface is trimmed. The trimming information is supplied in the 
global and local coordinate system of the surface to be trimmed. Eigure 1 shows an 
example of an IGES output. 

The basic idea of our approach is to use information provided in the IGES file 
directly for the definition of the geometry and other derived values such as the com¬ 
putation of the Jacobian and first and second derivatives. We are aware of only three 
papers that have addressed the problem of using trimmed CAD data for analysis. The 
first two papers ([1] and [2]) propose to use a regular grid of elements, that are defined 
by knot spans. A searching algorithm is employed, that allows to determine how the 
elements are transected by trimming curves. Eor the numerical integration trimmed 
elements are subdivided into triangular subregions. 


Available from http : / /diyhpl.us/ bryan/papers2/IGES5 — 3-forDownload.pdf 




Figure 1: Example of IGES output: Tunnel intersection and NURBS surface and trimming information for 
one tunnel 


The method is very general and can deal with extreme trimming cases such as 
multiple holes and cases where trimming curves are very close to each other. However, 
the implementation of the method is not trivial. The third paper [3] deals with the 
application of trimming to shell surfaces and with the specification of local loading. 
The method uses a reconstruction of knot spans and control points and is also generally 
applicable. 

The way our proposed method differs from published ones is that the implementa¬ 
tion in existing software is very simple. The application is restricted to cases where a 
unique mapping exists. In it’s current form, this excludes cases with multiple holes that 
have been analyzed in [1] but as will be shown, it is applicable to a number practical 
problems, especially in the main area of interest of the authors, tunneling. However, 
expansion of the method is possible. 

2. Trimmed NURBS, double mapping method 

For the explanation of the trimming method, we define two local coordinate sys¬ 
tems. One is the local coordinate system of the NURBS surface, used in standard 
isogeometric analysis, where the local coordinates v span from 0 to 1, the other is an 
additional coordinate system s^t, whose coordinates also span form 0 to 1. 
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Figure 2: Example of trimming. Left: Surface to be trimmed, Right: Two trimming curves with control 
points 


The new approach to trimming relies on the mapping of an area bordered by 2 
trimming curves and straight lines connecting the ends of the curves (defined in the t/, v 
system) into the s^t system. The global coordinates can now be expressed in terms of 
the local coordinates as: 


x{u,v) =x{u{s,t),v{s,t)) (1) 

This allows all operations, such as computing the derivatives of functions and the 
numerical integration, to be carried out in the sj coordinate system. A visual explana¬ 
tion of this procedure is provided in Figures 2 and 3. 



Figure 3: Explanation of mapping: right: s,t map, middle: u, v map, left: map 
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We start with mapping the two trimming curves defining the boundaries of the 
trimmed surface: 
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where RUs), Rn{s) are one-dimensional NURBS basis functions defining the trim¬ 
ming curves, are the number of control points and and the 

w,v coordinates of control points. The superscript refers to the curve number (//=top; 
/=bottom). As can be seen in Figure 3 these curves map as straight lines in the 5*,^ 
coordinate system. 

Next we use a linear interpolation between the curves to map the trimmed area: 


u{sj) = Ni{t) • u\s) +A^2(0 ’(3) 
v{s,t)=Ni{t)-v'{s)+N 2 {l)-v'\s) (4) 


with 

Ni{t) = l-t;N2{t)=t (5) 

Finally we map from the t/, v system to the global v,y,z system by 
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where v) are NURBS functions describing the untrimmed NURBS surface. 

A^B is are the number of control points in u^v directions and Xa^ij are control point 
coordinates. 

For the computation of the Jacobian matrix we need the first derivatives of x in 
terms of 

The first derivatives of x to 5* and t are given by 
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and 
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For Kirchhoff shell analysis we require the second derivatives also, which are given 
by; 
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The derivatives of u and v are given by: 
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The method can only be used in conjunction with the concept of geometry inde¬ 
pendent field approximation, which was first introduced in 14] and applied to practical 
problems in 15], 16], 17] 18] and 19]. 

This means that the description of the geometry is unchanged during refinement 
and that the field variables are approximated independently. For example, unknown 
displacements are approximated by 


u = 


b=la=l 


( 12 ) 
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where N^f{s,t) are basis functions defined in the s^t coordinate system, are 
parameter values and are the number of parameters in directions. The basis 

functions and the integration regions are defined in the local sj coordinate system and 
then mapped into the global system. Because the mapping affects the continuity of 
the functions in the global parameter space, care has to be taken in the choice of the 
functions and integration regions (see comments later). 


3. Implementation 

The method was implemented into isogeometric FEM and BEM software using 
Octave and the nurbs toolkit^. See [4] and [10] for more details about the program. 
The main mapping function is shown here: 


function [uv,vecu]= Map(Knotg,Coefs,tt) 


Maps points on trimmed surface from s,t to u,v 
Input: 

Knotg ... array containing knot vectors of all trimming 
curves 

Coefs ... array containing control points of all trimming 
curves 

tt (1,np),tt (2,np) ... vector with s,t-coordinates 
Output: 

uv(1,np),uv(2,np) ... u,v coordintes of points 
vecu ... Vectors in u,v directions 


for ncrv=l:2 

% extract knot values and control points of trimming curves 
[knotu,coefs]= Get_infoc(Knotg,Coefs); 

% mapping of trimming curves from s,t to u,v 
nurbs= nrbmak(coefs,knotu); dnurbs= nrbderiv(nurbs); 

[uv,vuv]= nrbdeval(nurbs,dnurbs,tt(1,:)); 
uvn(:,:,ncrv)= uv(:,:); vuvn(:,:,ncrv)=vuv(:,:); 
end 

% mapping of points on trimmed surface from s,t to u,v 
for np=l:columns(tt) 

Nl=l-tt(2,np); N2=tt(2,np); 

% compute u,v values 

uv(l,np)= Nl^uvn(1,np,1) + N2*uvn(1,np,2); 

uv(2,np)= Nl^uvn(2,np,1) + N2*uvn(2,np,2); 

% first derivatives 

vl(l)= Nl*vuvn(1,np,1) + N2*vuvn(1,np,2); 

vl(2)= Nl*vuvn(2,np,1) + N2*vuvn(2,np,2); vl(3)= 0; 

v2(l)= uvn(l,np,2) - uvn(1,np,1); 

v2(2)= uvn(2,np,2) - uvn(2,np,1); v2(3)= 0; 

vecu{1}(:,np)=vl; vecu{2}(:,np)=v2; 


^available for free download from http://octave.sourceforge.net/nurbs/overview.html 
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end 

endfunction; 


To implement the method insert the mapping function before the mapping to global 
coordinates. For example: 


[uv,vecu]= Map(Knotr,Coefstr,tt); 
[x,vect]= nrbdeval(nurbs,ders,uv); 



Symmetry plane 


Figure 4: Test example 1: Plate with a hole 


4. Test examples 

We have chosen two test examples, which are the same or similar to the ones used 
in the cited references on trimming. They relate to the FEM analysis of a plate and a 
Kirchhoff shell. 

4.1. Plate with a hole 

This is the exactly the same test example that has been used by Kim et al [1] to 
verify their trimming method. It relates to a plate with a hole. Exact solutions are 
available for an infinite plate with a hole. If the correct boundary conditions are applied, 
then a plate of finite extent can be analyzed and the solution compared with the exact 
solution. Figure 4 shows the geometry of the plate with the symmetry planes assumed. 
The size of the plate is 5x5 and the radius of the hole is 1. The elastic properties 
assumed are £" = 10^, v = 0.3. The plate is subjected to a tensile horizontal stress of 1 
at the right boundary. To simulate an infinite extent, the tractions computed from the 
exact solution were applied at the top. The trimming is explained in Figure 5. 
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Figure 5: Geometry definition of plate with a hole: Top left: untrimmed NURBS surface in global x,y system 
and Top right: trimming curves in local u,v system. Bottom: trace of basis functions along r = 1 in the local 
s^t and in the global x,y coordinate system 
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The definition of the 2 trimming curves is as follows: We start with a bilinear 
NURBS surface describing the plate without the hole and add two trimming curves. 
The plate is described by: 


Order: 1 

Knot vector: 0011 
Coefficients (x,y,z,weight): 
0 0 0 1 
0 10 1 
10 0 1 
110 1 


The trimming curves are defined by: 


% Curve 1: 

Order: 2 

Knot vector: 000111 
Coefficients: 

0 0.2 0 1 

0.2 0.2 0 0.707 

0.2 0 0 1 

% Curve 2: 

Order: 1 

Knot vector: 000.511 
Coef f dents : 

0 10 1 
110 1 
10 0 1 


Because one of the trimming curves has only a continuity, care has to be taken 
when choosing the basis functions for the description of the unknown before refinement 
and the regions used for the integration. In Figure 5 we show that in this case the 
integration region is split into two and the continuity of the basis functions is changed 
to match the one of the mapping function. In this Figure we also show the trace of 
the basis functions for the approximation of the unknown along t = \ in the local s^t 
coordinate system and in the global v, y coordinate system. 

We analyze three refinements and as expected, the convergence of the norm of 
the stress shown in Figure 6 is similar to the one reported in [1] for analyses without 
trimming. 

4.2. Kirchhojf shell 

This test example is similar to the one used by Schmidt et al in [3] to test their 
trimming method. It is the analysis of a Kirchhoff shell. Details of the implementation 
of isogeometric Kirchhoff shell analysis can be found in [11]. For the test we have 
implemented our mapping method into software developed as part of a Master Thesis 
[ 12 ]. 
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Figure 6: The 3 refinement stages investigated and convergence plot 
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Figure 7: Geometry definition of test case with 2 untrimmed NURBS patches, showing control points 


Two analyses are performed: one with two untrimmed and the other with two 
trimmed NURBS patches. 

Figure 7 shows the geometry definition with 2 untrimmed NURBS patches. The 
patches are fixed in three coordinate directions at the base. 

At the top only the translational degrees of freedom are connected (hinged con¬ 
nection). The shell is subjected to selfweight. Figure 8 shows the case where the two 
NURBS patches have been extended and then trimmed back to match the geometry of 
the untrimmed case. Both cases were analyzed and convergence achieved by order el¬ 
evation. The convergence of the maximum defiection for both cases is shown in Figure 
8 and as expected, no difference between the cases can be observed. 


5. Practical examples 

Here we show two practical examples that demonstrate that the proposed mapping 
method, despite its limitations, has practical applications. One example is in shell 
theory and is meant to demonstrate how easily the geometry can be changed. The 
second one relates to the boundary element analysis of a branched tunnel. 

5.7. Example 1: Arched Scordelis roof 

As a first example we show a variation on a classical example in shell theory, the 
Scordelis-Lo roof. The geometry of the problem is shown in Figure 9. 

The shell is subjected to self weight q and the properties are: 

• E= 4.32■ 10*^ kNjn? , v = 0 

• t = 0.25 , q = 90 kN/rtf' 

We change the geometry by introducing some arching and use three different trim¬ 
ming curves with increasing removal of material (Figure 9). 
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Figure 8: Two untrimmed NURBS surfaces with trimming curves, resulting trimmed surface and Conver¬ 
gence of maximum displacement 
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Figure 9: Scordelis-Lo roof: Problem definition and the trimming curves for the introduction of different 
degrees of arching (from left to right geometry 1 to 3) 
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Figure 10: Arched Scordelis roof: Convergence of the maximum displacement 




Figure 11: Contours of vertical displacement for geometries 2 and 3 
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It can be seen in Figure 10 that all cases converge. For moderately arched shells 
the maximum vertical displacement decreases but if more material is removed there is 
an increase. 

In Figure 11 we can see that for moderate material removal the maximum displace¬ 
ment is localized, which spreads over the shell as more material is removed. 

5.2. Example 2: Tunnel branch 

This example is meant to demonstrate the implementation of the mapping method 
in a Boundary Element code. The implementation of the isogeometric BEM is dis¬ 
cussed in [4] and of the mapping method in [10]. Here only results of the simulation 
are presented. More details can be found in [13]. 

Figure 12 shows the Rhino CAD model of a tunnel intersection and the geometrical 
information extracted from the IGES file produced by the program. This information is 
used directly for the simulation as follows: The trimmed surfaces are used to model 1/4 
of the problem as shown in Eigure 13. Two planes of symmetry (about the x-y and x-z 
planes) are assumed and specially developed infinite plane strain patches (see [4] and 
[10] for details) are used to simulate infinitely long tunnels. The tunnels are assumed 
to be excavated in one step in an elastic pre-stressed ground. 

Remark: Since the trimming curves are specified separately for each NURBS 
patch they will in general not match perfectly and there may be small gaps. The gaps 
can be reduced to a negligible size by specifying a small tolerance (option ’’units and 
tolerances” in Rhino). The other problem is that the parameter spaces of the trimming 
curves may not match at the interface. The implication of this for the BEM is that the 
collocation points, whose coordinates are computed using Gervilla abscissa in the lo¬ 
cal coordinate system of the trimmed NURBS, may not be in the same location. If this 
issue becomes significant the way of dealing with this would be by using discontinu¬ 
ous collocation (see for example [14]). Eor the example presented here the collocation 
point locations matched within 1 %, so continuous collocation was used. 

Using geometry independent field approximation, the basis functions for the def¬ 
inition of the unknown were refined until convergence was achieved. The displaced 
shape of the tunnel walls is depicted in Eigure 14. 

A comparison with a standard isoparametric BEM analysis is shown in Eigure 15 
and good agreement can be observed. 
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Figure 12: CAD model of tunnel intersection and extracted NURBS surface and trimming information 
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Figure 13: Geometry description for the simulation 



Figure 14: Deformed shape of the tunnel 



Figure 15: Comparison of vertical displacements along the intersection line 
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Figure 16: Example of a more complex tunneling situation 


6. Summary and Conclusions 

A novel approach to trimming was presented. The main advantage of the method 
lies in its simplicity and ease of implementation into existing software. The method is 
currently restricted to cases with two trimming curves. However, as has been demon¬ 
strated, practical applications exist, where this is the case. Our main interest lies in the 
application of the BEM to the simulation of underground excavations, where benefits 
of the method can be seen. The geometry of the underground excavation depicted in 
Figure 16 for example is ideally suited to the proposed algorithm. Using a conventional 
BEM simulation, it has been found difficult to generate a mesh from the CAD data, as 
elements with bad aspect ratios appeared and had to be repaired. Since the trimming 
information provided by the CAD program only involves two trimming curves that are 
not straight lines in this case, this application lends itself well to the proposed simula¬ 
tion approach. 

Extensions of the method are possible. In fact, Gordon-Coons patches (see [15]) 
provide a framework where four curves can be used to define a surface. One may 
thus construct an arbitrary patch with four trimming curves in the w, v-parameter space. 
Multiple holes may be considered by using a subdivision of the mapping area. 

Trimming plays an important role in realizing the dream of a direct connection 
between CAD and simulation, without the intermittent step of mesh generation. It is 
hoped that our contribution provides impetus for a - much needed - increase in research 
activity in this area. 
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